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We present a framework for efficiently performing Monte Carlo wave-function simulations 
in cavity QED with moving particles. It relies heavily on the object-oriented programming 
paradigm as realised in CH — h, and is extensible and applicable for simulating open interact- 
- ing quantum dynamics in general. The user is provided with a number of "elements" , eg 

pumped moving particles, pumped lossy cavity modes, and various interactions to compose 

O : 

£N1 . complex interacting systems, which contain several particles moving in electromagnetic fields 

of various configurations, and perform wave-function simulations on such systems. A number 
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of tools are provided to facilitate the implementation of new elements. 



INTRODUCTION 



u 



Based on our experience gained in recent years in Monte Carlo wave-function (MCWF) situ- 
ations of simple moving-particle cavity QED (CQED) systems performed with low-level codes 
, we have decided to summarise our know-how on the problem by developing a high-level 
framework for such simulations. The framework is highly modular and therefore easy to main- 
tain, relies solely on standard C++ programming techniques and therefore portable, and provides 
c3 an interface which is easy to use even for those not so familiar with the theoretical models of 

^ ! n 

Q" 1 moving-particle CQED ([5f| is a review of the theory involved). Meanwhile, thanks to the optimi- 



sation mechanisms of C++ compilers, we are safe to claim not to have noticeably lost in efficiency 
as compared to our previous low-level codes. Potentially, the framework is of good use for the 
quantum optics community. 

Simulating moving quantum particles presents many non-trivial numerical problems especially 
of stability l|, |4[. Hence, in the framework very careful numerics is needed. Accordingly, as 
discussed in App.|Aj we use a slightly modified version of the original MCWF algorithm (cf eg [(J) 
involving the use of adaptive step-size ODE steppers and interaction picture. 

At present the framework consists of three parts: The first part is the MCWF driver (Sec. Mil) . 
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which has only an abstract view on the open system to be simulated, represented by an abstract 
class. This abstract class stands at the origin of a class hierarchy consisting the second part of 
the framework (Sec. lIVp . Eg a system can be an element system or a composite system contain- 
ing several element systems. The aim of the hierarchy is to provide the user with tools to build 
composite systems from several elements, and to facilitate the implementation of such elements. 
Clearly, the first two parts stand quite independently of each other and are also generally appli- 
cable. As the third part of the framework several elements are provided at the lower levels of the 
hierarchy intended as building blocks for systems of moving-particle CQED (Sec. El). The building 
blocks are pumped moving particles, pumped lossy cavity modes, pumped two-level atoms, and 
interactions between them eg interaction between a cavity mode and a pumped particle moving 
along or orthogonal to it. This third part is independent of the first part, but not, of course, of 
the second part, the elements stemming from the same class hierarchy. 

For a given system on the highest level the user is required to write a simple driver program in 
C++ in which he/she defines the system to be simulated using the elements (selecting a number 
of free elements and interactions between them) and passes this system to the MCWF driver, 
which then evolves the system on a number of Monte Carlo trajectories. A description of the user 
interface and example drivers are given in Sec. HU 

Note that our approach here is quite different from the one presented in [?|]. That approach 
is built on a hierarchy of classes representing Hilbert space operators and state vectors, and the 
application of operators on vectors is defined. Operators acting on complex systems can then be 
built from elementary operators using direct product. A similar idea is implemented in the popular 
Quantum Optics Toolbox for MATLAB 8j. Consider for a moment how this approach could be 
applied for moving particles: In this case dealing with both operators x and p cannot be avoided. 
A moving-particle state vector can be stored in either representation (the state-vector object stores 
in which representation it is at the moment), and when the other operator is to be applied, an 
in-place Fast Fourier Transformation (FFT) is needed. However, as our experience shows, such a 
transformation always has numerical errors, which can disturb careful statistics. 

In our approach the user is provided with a much higher level interface, our classes representing 
whole physical systems instead of Hilbert space operators. This is certainly at the cost of flexibility, 
but our framework does not aim at such generality as the above, since it has been developed with a 
more concrete problem in mind, in particular, CQED with moving particles. For this given problem 
we consider our approach as more efficient than the above, since, as we will show in Sec. |V]we can 
completely avoid in-place FFT. 



In the following we first present the highest level of the framework, that is, the user interface, 
so that the reader can immediately get a feeling about our approach. Also, by reading Sec. HTlthe 
reader can in principle already use the framework, so this can be considered as a short write-up. 
This is followed by the long write-up, the presentation of the different parts of the framework. We 
include sections entitled "Desideratum" in which we indicate features that would logically belong 
to the given part, but are as yet missing because we have not yet needed them. These may easily 
be implemented in the future. 

Finally, in Sec. IVH we summarise our test runs performed with the framework. In the Appen- 
dices we describe our version of the MCWF method and the most important modules used in the 
framework. 

The source code contains more than 60 source files and a totality of about 4000 
lines, and is distributed in tgz format. It can be get either from SourceForge.net at 



http://sourceforge.net/projects/cppqed/ or directly from the authors. The framework has been 



tested under Debian GNU/Linux and RedHat Linux operating systems, in both cases the GNU 
C++ compiler has been used for compilation. 



II. THE USER INTERFACE 
A. Writing drivers 

The classes a user has to know about are listed in Tab. H together with the most important 
parameters, which will be explained further down in the text. The set of elements for systems in 
moving-particle CQED is explained in detail in Sec. El 

To ease the understanding of the framework's workings example drivers are given in Figs. [T] 
and [21 The driver in Fig. [1] simulates one single particle moving in a ring cavity, that is, two 
travelling-wave modes propagating in opposite directions. Both modes are lossy and one of them 
is pumped. In addition, the particle can also be pumped and scatter light from the pump into the 
modes. The driver in Fig. [5] describes two identical particles moving orthogonal to the axis of a 
single-mode cavity in a standing-wave pump field. 

The user has to choose an appropriate set of free systems and the interactions between them, 
and instantiate the corresponding Free and Interaction classes with the appropriate parameters. 
If two elements are exactly identical, only one object is needed. This is the case eg with several 
identical particles: one instant of the MovingParticle class stands for all of them (an example of 
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Ac, k, photonCutoff 


PumpedLossyMode 
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MovingParticle 
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PumpedMovingParticle 
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Interactions 




ParticleOrthogonalToCavity 


cavity, pumpedParticle, Uq 


ParticleAlongCavity 


cavity, (pumped)particle, £7o, r] e s, -Kcavity, cavityModeFunction 


ParticleCavity2D 


cavity, particle, pumpedParticle, Uq, r] e s, -Kcavity, cavityModeFunction 


ParticleTwoModes 


particleCavityl, particleCavity2 


IdenticalParticles 


(pumped)particle, iVp art i c ie, vector< partic i c )> 


Composite 


vector<SubsystemsInteraction> 


Subsystems Interact ion 


Inter act ion&, vector<subsystemSequentialNumber> 


HS_Vector 


dimension 


Trajectory 


\^${t = 0)), OpenSystemfe, seed, eps, dplimit 



Tabic I: Classes constituting the user interface of the framework, with the set of elements extendable in the 
future at will. Next to each class their most important parameters are listed, these are explained in the 
text. The Interactions take references to their subsystems as parameters — cavity is an instant of class 
(Pumped) LossyMode, (pumped) particle one of (Pumped)MovingParticle, while particleCavity one of 
MovingParticleCavity cf Sec. IV C II 



this can be seen in Fig. [2j). 

The Free objects are then to be (virtually) arranged into a sequence starting with number 0, 
and the user has to create a vector of Subsystemslnteraction class objects. The latter is a helper 
class for the Composite class, storing a reference to an Interaction and the sequential number 
of those Free objects between which the given interaction acts. Most interactions will be between 
two subsystems, but we have found cases with interactions between three or four subsystems (cf 
Sec. LY]). The IdenticalParticles class is an Interaction between all the particles, that is, an 
arbitrary number of subsystems in principle. 

When giving the sequential numbers the user has to remain consistent with the originally 
conceived sequence of the Free objects, and the order of the subsystems in an Interaction object 
is also important. Eg in Fig. [1] Line 19 instead of (pel, 1,0) it would be an error to write (pel, 1,2) 
because Free Nr. 2 is a PumpedLossyMode and not a MovingParticle, but also (pel ,0, 1) because 
ParticleAlongCavity is an interaction between a LossyMode and a MovingParticle and not vice 



versa. Such errors cause an exception during the construction of the Composite object. 

The free systems provide helper functions to prepare state vectors (of class HS_Vector) charac- 
teristic to the given system. Eg for a LossyMode object one can prepare a Fock state or a coherent 
state. This, together with the possibility of making up direct products of several state vectors, 
facilitates the preparation of initial conditions. Eg in Fig. Q] Line 26 we prepare a state in which 
the particle has a wave packet centred at position xO with momentum kO and spread xsig, the +K 
mode is in a coherent state with complex amplitude alpha, and the —K with beta. Here again, 
we have to comply with the our preconceived order of the Free objects in the sequence. 

As output such a program first summarises the parameters of the system, then at certain time 
instants (whose frequency is specified by the user) displays the time and the time step followed by 
a set of quantum averages specified in the element system classes. At specified time instants the 
whole state vector is displayed, but in practice this can be too big to store and gain information 
from. An example output is given in Fig. 

B. Desideratum 

With some effort the preparation of drivers could be made automatic, such that the user is 
presented with a higher level interface in which he/she specifies the system using some simple 
formal language, and then the framework writes and compiles the C++ driver corresponding to 
the system. A similar idea can be found implemented in the XMDS package ({J. 

III. EVOLUTION 
A. MCWF trajectories 

What we expect from a MCWF trajectory driver class (called Trajectory in our framework); 
what parameters does it need and what functionalities should it provide? 

First we need to represent the state vector of the system. The most straightforward representa- 
tion is a complex packed array (CPA), that is, a real array, in which the real and imaginary parts 
of the state- vector amplitudes are placed in alternate neighbouring elements. In our framework the 
low-level notion of a CPA is furnished with an interface class called HS_Vector (for Hilbert-space 
vector) supplying the operations we expect for a vector of a Hilbert space. These include algebraic 
(vector-space) operations including direct product of several vector spaces, metric operations, and 
both low and high level access to amplitudes. When instantiating a Trajectory the initial condi- 



tion of the system has to be given in an appropriate instant of HS_Vector and this is eventually 
replaced by the driver when evolving the system. 

Every system must supply an interface towards the trajectory driver containing the operations 
needed to perform a MCWF step on the system as described in App. S This interface is the ab- 
stract view the driver has on the system to be simulated. In the present framework such an abstract 
system is represented by an abstract class called OpenSystem. The hierarchical implementation of 
this interface for more and more concrete systems constitutes the main part of the work presented 
here and is described in Sec. IIVI In the C++ implementation of the object-oriented paradigm, an 
abstract class cannot be instantiated but can be referred to by a reference (a pointer), to preserve 
run-time polymorphism. Hence, a Trajectory object takes a reference to an OpenSystem. 

An ODE integrator and a random-number generator are needed to perform Step 1 and 2 of 
an MCWF step, respectively. These are also wrapped into interface classes called Evolved and 
Randomized, respectively. At the moment, these classes are implemented using the Gnu Scientific 



Library (GSL) lCj], but here a user is free to choose his/her own favourite library (eg Numerical 
Recipes) or even hand-crafted code. To better localise object creation, only "factory" objects for 
these classes are passed to the Trajectory object (for a description of the factory-class and other 
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programming techniques appearing in this paper see 

Other important parameters are the highest allowed jump probability dplimit and the relative 
precision for the ODE stepper eps. 

The class supplies a member function called Step to perform one adaptive-stepsize MCWF step 
on the system as follows: 



1. Invokes the ODE stepper to evolve the state vector according to Eq. (|A2j) for a suitable time 
interval dtdid. H n n for the system is taken from the OpenSystem class. 

2. Performs the additional (exact) part of the evolution as \^>(t + 6t)) = U^ 1 (5t) \^i(t + 6t)). 

3. Examines whether a jump should be made. For this it uses a random number, dtdid, and 
a system-specific jump function again taken from the OpenSystem class. 

4. The ODE stepper supplies a time step dttry which is likely to work for the next step. 
The driver examines whether the jump probability would have overshoot dplimit were it 
calculated with dttry instead of dtdid. If this is the case, dttry is reduced. 

5. Calculates and communicates towards the user physical properties of the system at the given 
time instant, such as the state vector itself and/or important quantum averages — exactly 



what is again taken from OpenSystem. 

A number of helper functions are provided to take not only a step but evolve a whole trajectory 
or an ensemble average of trajectories. 

B. Desideratum 

Other methods of wave-function simulation of open systems can be straightforwardly added to 
the framework, although the OpenSystem interface may need to be extended. These include the 
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quantum state diffusion method [12|, and the orthogonal quantum jump method [13]. It would be 
advisable to keep a common interface for the different drivers, so that the same helper functions 
work for all of them. 

Wave function simulations can very efficiently be done parallel. With additional helper functions 
parallel execution can be easily implemented. 

IV. SYSTEM HIERARCHY 

Every class derived from OpenSystem is an OpenSystem, features the same interface, and hence 
can be passed to the trajectory driver. 

As indicated in Fig. 0] an OpenSystem is either Composite or Element system. Element systems 
can be used as building blocks to compose composite systems. One may wonder why derive also 
Element from OpenSystem when elements are simple systems with known behaviour, so that one 
is unlikely to wish to simulate such systems. The answer is that one may wish to simulate them 
for testing purposes when implementing a new Element class. Also, this way quite an amount of 
code can be reused. 

An Element, in turn, can be either a Free system or an Interaction of such systems. We 
emphasise the fact that an Interaction is also an Element, and hence an OpenSystem. One is 
even less likely to wish to simulate only the interaction part of the dynamics without the free 
systems: The reason for this arrangement is again code reuse. 

We note that we had considered the alternative design depicted in Fig. [SJ Here, there is a 
very clear distinction between system that use interaction picture and those that do not. In 
many sense this design is more logical and attractive, since it grasps better the structure of the 
problem. However, it involves the use virtual bases, consisting a slight efficiency overhead, and, 



more importantly, a bigger overhead in the complexity of the code. We therefore eventually resorted 
to the first simpler design for the testing phase. 

The design we have found ultimately useful is, however, the one depicted in Fig. El This one 
unites the advantages of the previous two, without the overhead of virtual bases. This design is 
uncompromising in the sense that it is very clearly expressed which virtual functions a class at the 
lower levels of the hierarchy has to implement. 

Although the underlying design in our framework is this last one, in the following, for the 
sake of simplicity, and to ease the understanding for those not so familiar with object-oriented 
programming, we go on presenting the framework as if the underlying design was the first one. 
The differences are purely technical throughout. 

A. OpenSystem 

The OpenSystem class is not a purely abstract one, since it has one data member: the dimension 
of the system — a parameter every quantum system has in common. In addition it features 
a number of virtual functions (function prototypes) which enable the driver class to perform a 
MCWF step as described in Sec. IIII Ai Eg the (non-Hermitian) Hamiltonian of the system is 
implemented by the function 

void H ( double t, const double* Psi, double* dPsidt, const CPA_View& V ); 

The first three arguments are the expected ones: time, an array for the state vector \fy), and one for 
the state- vector derivative d \^f) jdt. It is the last parameter that needs some explanation. Since an 
OpenSystem can be an Element system, it must be prepared to be embedded into a complex system 
as a subsystem. If so, to be able to perform the operation on the state vector of the whole system, 
H must have some information about the embedding complex system. As explained in App.[Bl this 
information can be condensed into a set of array slices, which set, in turn, is implemented by a 
class called CPA_View in our framework. 

The other important virtual member functions are U, J, and Display, which take care of Phases 
2, 3, and 5 of a MCWF step as described in Sec. IIII Al respectively. They all take arguments one 
would expect them to, plus a CPA_View. 

A further important virtual member function is called HighestFrequency, and returns the 
highest characteristic frequency in the system's time evolution — a measure what every dynamical 
system is expected to have. This is needed by the Trajectory driver to initiate the ODE stepper: 



adaptive step-size ODE steppers need a good guess for the initial time step to try, which is derived 
by the driver from the highest characteristic frequency of the system. 

B. Element 

At the level of DpenSystem the functions H, U, J, and Display are virtual functions because we 
can not tell what these functions are to do for a general OpenSystem. 

An Element system will be mostly embedded into a complex system as a subsystem. As 
explained in detail in App.[B]to calculate eg the Hamiltonian it has to iterate over the state- vector 
slices contained by its CPA_View, which corresponds to iterate over all the possible combinations 
of the quantum numbers of other subsystems — the "dummy" quantum numbers from the given 
subsystem's point of view, and call the same function on the corresponding slice. Function H is 
implemented accordingly, and class Element hence features the virtual function 

void H_elem ( const double* Psi, double* dPsidt, const CPA_Slice& S ) const; 

Note that the time argument is not passed over to H_elem. The time dependence of the original 
Hamiltonian H is rather taken care of by another virtual function H_update, which updates the inner 
state of the object if it does not correspond to the given time instant. With this method much 
calculation can be saved when the same object is used to describe several identical subsystems. 

Note that Element is also an abstract class because although it implements function H from 
OpenSystem, it declares new virtual functions, which must be implemented further down in the 
hierarchy. 

J and Display are implemented along similar lines as H, in both cases new virtual functions are 
declared. Eg for J we need a function J_dpoverdt which calculates the probability of a jump per 
unit time in the given subsystem, and a function J_elem which actually performs the jump on a 
given state-vector slice if required. 

U is not implemented by Element. An Element can be Free or Interaction. U represents the 
part of the dynamics which can be exactly solved, that is, the part of the Hamiltonian which can 
be diagonalised. This is possible for some free systems, but not for interactions. Therefore U is 
implemented only in class Free, along exactly the same lines as H in class Element. 

Interactions may affect the parameters of frees. A straightforward example for this is a cav- 
ity mode whose resonance frequency is shifted when interacting with an atom. Hence, class 
Interaction features a virtual function called FreesAdjust, which performs the required mod- 



ification in the parameters of the free systems. It is important to note that this is done at the 
construction of Composite rather than at the construction if the given Interaction. Indeed, at 
the construction of the interaction we do not yet know how many times it will be applied: this 
becomes clear only when we already know the layout of the whole composite system — in the 
above example the cavity frequency has to be shifted twice if there are two atoms instead of one. 

Not every element has to implement all the virtual functions declared in class Element. Eg we 
can easily imagine free systems whose dynamics can be exactly solved. In this case the coherent 
evolution is completely taken care of by U, hence H_elem and H_update need not be implemented. 
An other common case is when an element's dynamics is purely coherent. In this case the func- 
tions connected to J are not implemented. An interesting case is that of IdenticalParticles cf 
Sec. IV C 4( which can be considered the extreme: this class exists solely to perform calculations in 
occupation-number representation, and implements solely the functions related to Display. 

C. Composite 

A very important task of class Composite is to keep track of its elements (frees and interactions) 
and their CPA_Views. The calculation of the CPA_Views takes place already at the construction of 
the Composite object. 

Composite is a concrete type, so that it has to implement all the virtual functions of its parent 
class DpenSystem. Eg H is implemented as calling successively the H of each element with the 
CPA_View of the given element. For this to work, it is important that the H functions of the 
elements add their contribution to dPsidt rather than replace it. Hence with the successive calls 
the contributions of elements add up, according to the model (|B1D . 

The implementation of the composite U and Display is rather similar, only J needs a bit more 
elaboration, since here the element Js should not be performed one after the other, but a choice 
has to be made as to which one (if any) to perform. The interested reader should refer to the code 
to see how this is implemented. 

D. Desideratum 

It is an interesting possibility, and one whose implementation should not be too difficult in the 
framework to allow composite systems to be elements of even more composite systems. This would 
be useful eg to facilitate the simulation of several atoms of complex structure. 



V. EXAMPLE: MOVING PARTICLES IN CAVITY 



A. Theory 

Let us consider a single pumped two-level atom interacting with a single pumped lossy cavity 
mode. Our units are chosen such that H = 1. Using the Jaynes-Cummings model to describe the 
arising interactions, the Hamiltonian for such a system reads (a is the cavity field operator, the as 
are the atomic internal operators, r and p are the atomic position and momentum operators) 

2 

H = -Ac at a+i (rjat - rj*aj + ^- - A A <J Z + i (r^(r)a - r?t(r)cr^ -i (^g(r)a^a - g*(r)ataj , (la) 

where the terms describe free field, pumping of the mode, atomic external and internal degrees of 
freedom (free and pumped), and atom-mode interaction, respectively. The Liouvillean reads 

C P = k (2apa) - [a+o, p] J + 7 ^2 J d 2 u N(u) ae~ ikAUr p e ifcAU V - [aV, p] + J , (lb) 

where the first term describes cavity decay and the second one atomic spontaneous emission. The 
second term contains momentum recoil due to spontaneous emissions. The unit vector u is the 
direction of the spontaneously emitted photon, and iV(u) the direction distribution characteristic 
to the given atomic transition. 

The operator (jlbj) conforming with Eq. (jAlh we can immediately read the necessary jump 
operators for this system. There is one for cavity decay and an infinite set parametrised by u for 
atomic decay: 

J c = V2kq, J a (u) = ^ie- iKur a. (2) 

We introduce Zq = k — iAc, Z\ = 7 — iAj^. The non-Hermitian Hamiltonian is obtained by 
replacing Ac with iZq and Aa with iZ\ in Eq. (|lal) . 

In the limit of large atomic detuning Aa the atomic internal degree of freedom a can be 
adiabatically eliminated, as described in Refs. 

_ g(r)a + r] t (r) 

~ ^Aa — 7 1 ' 

In this limit the atomic spontaneous emission can be neglected in most cases of interest. We will 
resort to this approximation to simplify the discussion. Putting 7 = leaves us with only one 
jump operator 
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Tabic II: Summary of the free elements' functionality, fully exposed in the text. f|- indicates that the given 
function is inherited from the parent class. TV = a^a is the photon number of the mode. 



We plug ([3]) into (fTaj) . We take <?(r) = gf(r) and r/ t (r) = 7?t,C(r), and assume that g and r]t are real 
(the possibility of their being complex is investigated in l4j). We obtain the following effective 
non-Hermitian Hamiltonian 



H, 



P 



+ sign(C7 )V^o%ff /*(r) C(r) o+ + h.c. , (4b) 



eflf 



with U = \g\ 2 /A A , i] c fi = |r/ t | 2 /A A . 

The following set of elements realizes the system An important restriction whose reason 
will become apparent later in this section is that the mode functions are restricted to one dimension 
and either standing- or travelling- wave modes: 



/(r), C(r) = m(£) = { 



sin(i^) 
cos(K£) , £ = x,y,z. 



(5) 



B. Free elements 



These classes implement H_elem, H_update, J_dpoverdt, J_elem, and the functions connected 
to Display: Average and AverageProc from parent class Element, and U_elem, U_update from 
parent class Free. Their functionality is summarised in Tab. HH 



1. (Pumped) LossyMode 



Class LossyMode implements the dynamics of a free lossy (cavity) mode. Its parameters are 
the detuning between the driving and the cavity resonance Ac, the cavity decay rate k, and the 
photon number cutoff. 



The non-Hermitian Hamiltonian can be diagonalised exactly, so that H_elem needs not be im- 
plemented while U_elem is implemented as applying U(t) = exp (— Zcta^a) on the state vector 
slice. 

A PumpedLossyMode has the additional parameter r]. Here only H_update and H_elem needs to 
be implemented to apply the Hamiltonian 

Hi(t) = iU' 1 ^) (not - v*a) U(t) = i Ue Zct (J - r]* e - Zct aj . (6) 

Since pumping does not affect the remaining part of the dynamics, all the other functions are 
exactly the same as for LossyMode, and PumpedLossyMode indeed has access to these functions: 
"inherits" them from the parent class LossyMode. This is the reason why in the class inheritance 
hierarchy in Fig.[3]PumpedLossyMode stems from LossyMode. Clearly, this technique can be applied 
to reuse a lot of code, and has indeed been applied throughout in our framework. 

2. (Pumped)MovingParticle 

A similar relationship exists between MovingParticle and PumpedMovingParticle. 
MovingParticle implements the dynamics of a free quantum mechanical particle moving in ID, 
with Hamiltonian H = p 2 j(2\i). This Hamiltonian is most conveniently implemented in momen- 
tum basis. For the numerics the momentum basis must be discrete, which amounts to some finite 
quantisation volume (length). Our choice of units is such that the smallest momentum is Ak = 1, 
that is, the quantisation length is 2ir. It is easy to see that the use of discrete momentum basis 
entails periodic boundary condition at the borders of the quantisation length. The parameters 
are the recoil frequency lu tcc = % Ak 2 j (2[i) = 1/(2//) and the spatial resolution. The latter has to 
be an integer power of 2 to be able to perform radix-2 FFT on the state vector. With our units 
H = tv rcc k 2 with operator k = p/(hAk) = p. 

The Hamiltonian is diagonal in the momentum basis, and is quadratic in the momentum. 
According to our experience, the second property makes it essential to use interaction picture 
because the quadratic growth of the frequency is too quick for the stepper routine and results in 
instabilities for the higher momentum components. 

According to this discussion class MovingParticle implements U(t) = exp (— iuj vcc t A; 2 ) , which 
is diagonal in momentum basis. The quantum averages calculated and communicated towards 
the user are: (fe), (/c 2 ) — (A;) 2 (proportional to the kinetic temperature of the particle), (x), 
\l (x 2 ) — (x) 2 . This means that at each call of Display for the class, a Fourier transformation 
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Table III: Summary of the interaction elements' functionality. N = a^a is again the photon number, 
A = sign(J7 )\/C / o??cff- 

has to be performed on a copy of the state vector to calculate the averages of operator x. This is 
done using the radix-2 FFT routine supplied by GSL, but here again the user is free to use his/her 
own favourite routine. We emphasise, however, that the time evolution is performed purely in mo- 
mentum representation, nor is our Trajectory driver prepared to perform FFT during evolution. 
When FFT is performed at all, it is on a copy of the state vector, not an in-place transforma- 
tion. Hence, we avoid numerical errors accumulating in the state vector, and also save the inverse 
transformation (although we lose time by copying). 

PumpedMovingParticle implements the Hamiltonian Hj(t) = r) c fiU~ l (t) |C( r )| 2 U(t). This has 
to be done in momentum space as well, therefore it pays to choose £(r) such that it be easy to 
calculate its action on the state vector in momentum space. This brings us back to the restriction 
(jSJ): the action of e lK ^ is very easy to calculate as it simply amounts to a shift by K in momentum 
space. For £(r) = e ±lK ^ the Hamiltonian is constant, while for £(r) = sin(K^), cos(K^), it is 
proportional to =F cos(2K^)/2, respectively, after dropping the constant term. This gives 

Hl{ t) = T M (t) cos{2 K0 U(t) = T 1] f (e-^ec(tf-fc) e 2i^ + e -w^ c {K+k) e -2iK^ (7) 

It becomes clear how huge we gain by using interaction picture in this case. The Hamiltonian 
is time dependent now, but the oscillation frequency grows only linearly with k instead of the 
quadratic growth mentioned above. 

C. Interaction elements 

The functionality of these classes is summarised in Tab. IIIIl 



1. Particle(Orthogonal/Along)Cavity 



These classes implement the interaction Hamiltonians between a cavity mode and a particle 
moving in ID, either in a direction orthogonal to the cavity axis, or the direction along it, respec- 
tively. 

Hence, ParticleOrthogonalToCavity implements 



which describes atomic stimulated absorption of a photon from the atomic pump and stimulated 
reemission into the cavity mode or vice versa. ParticleAlongCavity implements 



where the first term describes atomic stimulated absorption from the cavity mode followed by 
stimulated reemission into the same mode. 

These Hamiltonians are also implemented in interaction picture. Note that the first Hamiltonian 
is formally identical to the second term of the second Hamiltonian. Therefore it pays to implement 
this term already in a higher level in the hierarchy, so that both of these classes have access to it. 
This is done by the class MovingParticleCavity which, as we see in Fig. 0] is a parent class of 
both. 

ParticleAlongCavity is either instantiated with a MovingParticle and an explicitly supplied 
parameter etaef f , or with a PumpedMovingParticle, in which case the etaef f parameter is taken 
from this latter class. The first case describes the situation when the particle pump is aligned 
orthogonally to the cavity axis, while the second case when it is along the axis, so that the particle, 
which is also moving along the axis, feels the pump potential as well. 

The virtual function FreesAdjust defined in Interaction is implemented so that the cavity 
frequency is shifted by the interaction with the particle. In the orthogonal case this is fairly straight- 
forward: the shift Ac — ► Ac — Uq is applied. In fact, the user has the choice whether it should be 
applied or not, in the latter case Ac stands for the shifted frequency. With ParticleAlongCavity, 
the situation is somewhat more involved because the shift depends on the cavity mode function: 
for /(£) = e ±lK ^ the shift has to be done by Uq, while in the /(£) = sin(K£), cos(K^) case by 
Uq/2 since in this case the first term of the Hamiltonian (|8b|) reads Uq/2 (1 =p cos(21f£)) a^a. 




(8a) 




(8b) 



2. ParticleCavity2D 



This class is an Interaction between three subsystems, and implements the Hamiltonian 

H = U |/(6)| 2 « f a + siga(U )^U^ (/*(&) C(6) a) + /(&) C*(6) «) , (9) 

which describes the situation when the pumped particle is moving in two dimensions. One of the di- 
mensions is taken care of by a MovingParticle class and the other one by a PumpedMovingParticle 
class — as mentioned above these classes implement one single spatial degree of freedom each. 

3. ParticleTwoModes 

It is easy to see that if we have several cavity modes then instead of ([3]) we have 

a oc ^0j(r)di +r/ t (r). (10) 

i 

In the effective Hamiltonian (|4b|) this creates terms like 

Focr(6)/(6)4«2+ii-c., (ii) 

which describes atomic stimulated absorption of a photon from one mode and stimulated reemission 
into the other mode. 

This cannot be described with the classes we have so far, so we need one more class 
ParticleTwoModes to cover this case as well. This closes our set of classes needed to build com- 
posite systems of an arbitrary number of (pumped) moving particles and (pumped) lossy cavity 
modes of different spatial configurations complying with the model @. 

ParticleTwoModes is an interaction between four subsystems, but the two spatial degrees of 
freedom can be the same. This describes the case of a linear cavity sustaining two modes and one 
particle moving along it. 

4- IdenticalParticles 

An interesting feature of our framework is that if we have several identical particles, it is 
very easy to switch between their being considered as bosons or fermions, or even distinguishable 
particles. All we have to do is to prepare the initial condition with the appropriate symmetry with 
respect to the swapping of two particles. This symmetry is then conserved during evolution. 



If we consider our particles as indistinguishable, we might want to perform calculations in 
some occupation- number basis. This is facilitated by the IdenticalParticles class, which is an 
Interaction between several identical particles, which are therefore described by one single object 
of class (Pumped) MovingParticle. At its construction, an IdenticalParticles takes a reference 
to such a particle object, the number of particles, and a set of single-particle state vectors. It 
then constructs the occupation number basis and Display is implemented such that the complex 
amplitudes in this basis are calculated and communicated towards the user. Of course this makes 
sense only if the single-particle state vectors are pairwise orthogonal. 

Eg for two particles and two state vectors \<fii) and |<fo) the occupation-number basis for bosons 
looks like 

|2,0) = |^i)<g>|&>, (12a) 
|1,1) = (\4>i) ® \(h) + 102) ® \4>i)) , (12b) 
|O,2) = |02>®|02>, (12c) 

and Display then displays the complex amplitudes (2,0 | X I / ), (1, 1 1 and (0,2 | 

In the case of indistinguishable particles it makes no sense to calculate the quantum aver- 
ages for each of them separately because due to the symmetry all will be equal. Therefore, 
IdenticalParticles implements FreesAdjust such that the Display of the particles is switched 
off, and taken over by IdenticalParticles. 



D. Desideratum 



We note that the above description of IdenticalParticles reflects the "ideal state" of the 
class, which allows it to be used completely generally. Clearly, for several particles and single- 
particle states the implementation of this involves an amount of combinatorics, and has not yet 
been done. Instead, in the first release of the framework IdenticalParticles is an interaction 
between two atoms, and calculates {n\n2}, where n\ is the number of particles at x < and 
n2 at x > 0. Why this is useful in some cases is explained in Q]. Of course, this restriction of 
IdenticalParticles does not mean that the framework can not be used to simulate as many 
particles as wanted. 

Atomic spontaneous emission is not implemented. In the above discussed model, where the 
atomic internal dynamics is eliminated, the implementation of this is rather involved, eg the jump 
operators have to be implemented by the interaction classes since they contain both operators x 



and a A physical problem with the spontaneous emission is that in the far detuned regime its 
rate is given by Tq = 7 g 2 / A\, which is much smaller than the other frequencies of the system. It 
therefore adds a new, very slow relaxation time scale to the system, which makes the simulations 
very long, practically unmanageable. 

It is interesting to note that when implementing spontaneous emission, class 
IdenticalParticles gains physical significance: it has to ensure that the particle jump 
operators do not modify the state vector's symmetry with respect to particle exchange. 

The next step in the development will be the addition of the two-level atom to the framework. 
This entails a number of new interaction elements, eg the term i (r]*(r)a — ( r )c^) of Hamiltonian 
(|la|) will be an Interaction between a two- level atom and one or several spatial degrees of freedom 
(MovingParticles). 



VI. TEST RUNS 



Testing is difficult in our case because the behaviour of the system we aim to simulate, that 
is, the coupled open quantum dynamics of several particles and lossy cavity-field modes is largely 
unknown, and constitutes an extremely rich area of active physical research — the framework is 
intended tool for this research. 

Of course, utilities like HS_Vector, Evoled, Randomized, and maybe even Trajectory can be 
tested separately. Free elements should not present too much problem either. Interactions are, 
however, more problematic. 

Our principle for testing interaction elements was to find parameter regimes where the action 
of one subsystem on the other(s) is very strong, but the back-action is negligible. 

As an example, imagine a very massive pumped particle moving quickly in a direction orthogonal 
to a cavity. The particle is initially prepared as a very well localised wave packet. The pump is 
weak, so that the atom does not feel any potential, but the coupling to the cavity mode is strong, 
although not strong enough to create a big field that would act back on the atom. In this case 
the cavity field is weak, but is very sensitive to the position of the atom, on the other hand, the 
atom does not feel the field at all. If the particle is quick enough, it can travel several pump 
wavelengths before its wave packet spreads noticeably. The cavity decay rate k is big enough so 
that the field follows adiabatically even this quick atomic motion. In this case in the initial phase 
of the dynamics the cavity field is almost a classical field scattered by an almost classical point-like 



particle. This field we can calculate explicitly: 



(a) = (13) 

— Uq+IK 

where C, is the pump mode function, and x is the position of the atom. An example for such a test 
run is displayed in Fig. [7J 
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Appendix A: DESCRIPTION OF THE MCWF METHOD 



The MCWF method 



aims at the simulation of open quantum systems based 



on a stochastic ("Monte Carlo") state vector. In terms of dimensionality, this is certainly a huge 
advantage as compared to solving the Master equation directly. On the other hand, stochasticity 
requires us to run many trajectories, but the method provides an optimal sampling of the ensemble 
density operator so that the relative error is inversely proportional to the number of trajectories. 

The optimal sampling is achieved by evolving the state vector in two steps, one deterministic 
and one stochastic (quantum jump). Suppose that the Master equation of the system is of the 
form 

9 = \ [P ' F] + Cp ~ \ [P ' H ^ + ^ ( Jm/ ° J ™ " 2 i J ™ Jm ' P \ + ) ' (A1) 

the usual form in quantum optics. At time t the system is in a state with normalised state vector 
|V?(t)). To obtain the state vector at time t + 5t up to first order in 5t: 

1. The state vector is evolved according to the non-unitary dynamics 

ih^ = H nU \V) (A2) 

with the non-Hermitian Hamiltonian 

H nU = H-^^2j} n J m (A3) 



to obtain (up to first order in Si) 



Since -ff n H is non-Hermitian, this new state vector is not normalised. The square of its norm 
reads 

+ St) I *nH(* + St)) = ( 1 + ( 1 " Mt)) =1-Sp, (A5) 



where 5p reads 

8p = 8t l - (tf (t)| tf nH - |*(t)) = s Pm, (A6a) 

m 

5p m = ft J^J m |¥(t)> > 0. (A6b) 

Note that the time step 5t should be small enough so that this first-order calculation be 
valid. In particular, we require that 

6p<l. (A7) 
2. A possible quantum jump with total probability 5p. For the physical interpretation of such 



a jump see eg Refs. [a, ll7|. We choose a random number e between and 1, and if 5p < e, 
which should mostly be the case, no jump occurs and for the new normalised state vector at 
t + St we take 

Mt + St))=^f\ St ». (A8) 
Vl - op 

If e < 5p, on the other hand, a quantum jump occurs, and the new normalised state vec- 
tor is chosen from among the different state vectors J m \ ^(t)) according to the probability 
distribution IT m = 5p m /Sp: 

mt + 5t))=Vsi Jm ]^- )) . (A9) 

Obviously, however, we can and must do much better than this. Indeed, assume that for some 
time no quantum jump occurs, and we perform Step 1 several times consecutively. This would be 
equivalent to evolving the Schrodinger equation with the most naive first order (Euler) method, 
which is known to be unstable and hence fail in most cases of interest. In our framework, we choose 
to use instead an adaptive step-size ODE routine, usually the embedded Runge-Kutta Cash-Karp 
algorithm In this case the time step is intrinsically bounded by a precision requirement in 

the ODE stepper, but also by the condition (1A7I) . which is taken care of by our MCWF stepper. 
Since in the ODE we are now much better than 0(St), the renormalisation of the state vector is 
performed exactly rather than to 0(St) as in Eq. (|A8D . 



In many situations it pays to use some sort of interaction picture, which means that instead of 
Eq. (IA2h we strive to solve 

ih ^ = U ~ 1 { H ^ U - ih lf)^ l) > (A10) 

where = U~ l Note that U can be non-unitary. The two pictures are accorded after 
each time step, ie before the time step = (^(t)) and after the time step the transformation 

\*5f (t + 5t)) = U(5t) \*5>i(t + 5t)) is performed. This we do on one hand for convenience and for 
compatibility with the case when no interaction picture is used, but on the other hand also because 
U (t) is non-unitary and hence for t — > oo some of its elements will become very large, while others 
very small, possibly resulting in numerical problems. It is in fact advisable to avoid evaluating 
U(t) with very large t arguments. 



Appendix B: INTERACTING SYSTEMS STATE VECTOR SLICES 

The main objective of the development of the present framework was to allow users to com- 
pose composite systems at will from elementary systems and interactions already provided in the 
framework, and perform simulations for these composite systems. We can think of quantum op- 
tics: several atoms of different structure interacting with light fields or cavity modes. A concrete 
example is given in Sec. [Vj 

Let us consider what we expect from an element of such a composite system. This element will 
be a class, containing all the necessary parameters specific to the given elementary system, and 
featuring eg a function which calculates the effect of the free elementary-system Hamiltonian H a t 
on a state vector. The Hamiltonian H for a composite system of N subsystems in terms of this 
Hamiltonian reads 

H = H + • • • + H at + ■ ■ ■ + H N + F interaction , (Bl) 

The action of the elementary Hamiltonian H at on a state vector \^) expanded in a basis specified 
by some quantum numbers {i n } n= o n can be written as 

({Un=0... J vl^at|^) = EKt m ). . fu M, '.Y * • (B2) 

Jat 

Since at the time of developing the class of the given elementary system we do not know in which 
environment it will be embedded, we expect the very same piece of code to work independently of 
the environment. On the other hand, it has to know something about the environment because as 



we see in Eq. (|B2p the multiplication by the matrix of H^ em has to be performed for all possible 
combinations of the "dummy" quantum numbers {«n} n ^ a r 

The state vector is ultimately stored as a one dimensional array (a CPA) no matter how complex 
the system is, and the quantum numbers {i n } n= Q n are ma PP e d to a one dimensional index by 
the indexing function 

N N 

I(i , . . . , i N ) = } \i n \\ d m , (B3) 

n =0 n+1 

where d denotes the dimension of the subsystem. Hence, the information needed by H at about 
the environment can be condensed into the concept of array slices, which, in our framework is 
implemented by the CPA_View class. For a free system, a CPA_View class consists of an array 
firsts which contains the indices I(io, ■ ■ ■ , i&t = 0, . . . , in) for all the possible combinations of the 
dummies {i n } nj L at and an integer stride = flat+i dm- 

To each element of the array firsts of a CPA_View there corresponds a CPA_Slice which contains 
one single index first and the integer stride. One can say that CPA_Slice is the iterator type 
of CPA_View. The index corresponding to a subsystem quantum number i a t for a given set of the 
dummy quantum numbers can then be calculated from the slice alone as 

I fiat! {^n} n ^ at ) = first + stride x i at . (B4) 

All the environment-independent implementation of H at and eventually that of every operator 
acting on a subsystem at of a composite system has to see from the environment is a CPA_View. 
Having received a CPA_View as a parameter all an elementary Hamiltonian H at has to do is to 
iterate over the dummy indices condensed into firsts and apply the same matrix i?at cm on t ne 
state-vector slice specified by the corresponding CPA_Slice. This concept is realized by H and 
H_elem, cf Sec. ITVBl 

CPA_View is essentially an array of CPA_Slices, we just save resources by storing stride, which 
is the characteristic of the given subsystem embedded in the given environment, only once. 

As discussed in Sec. IIVB[ interactions are also "elements" in our framework. An interac- 
tion Hamiltonian operates on several subsystems, therefore its CPA_View has to contain as many 
strides, each corresponding to the stride characteristic for the given subsystem in the given em- 
bedding environment. 

The concept of array slices, the relationship between CPA_Slice and CPA_View, and the fact 



that a CPA_View represents a way of looking on the state vector is further exposed in Fig. [HJ 
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#include "PumpedLossyMode . H" 

1 #include "ParticleAlongCavity .H" 

2 #include "ParticleTwoModes . H" 

3 #include "Composite . H" 

4 #include "Trajectory .H" 
5 

6 int main(int, char*) { 
7 

8 // Instantiate Frees 

9 MovingParticle p(omrec,f in) ; // FreeO 

10 LossyMode mPlus (DeltaC, kappa, cutof f 1) ; // Freel 

11 PumpedLossyMode mMinus (DeltaC , kappa, eta, cutof f 2) ; // Free2 
12 

13 // Instantiate Interactions 

14 ParticleAlongCavity pcl(&mPlus,&p,UO,etaeff ,K, Plus) ; // exp(iKx) mode (Plus) 

15 ParticleAlongCavity pc2(&mMinus,&p,U0,etaef f ,K, Minus) ; // exp(-iKx) mode (Minus) 

16 ParticleTwoModes ptm(&pcl ,&pc2) ; 
17 

18 // Instantiate Composite 

19 vector<Subsy st ems Inter act ion> i (1 , SubsystemsInteraction(pcl ,1,0)); 

20 i.push_back(SubsystemsInteraction(pc2,2,0)) ; 

21 vector<size_t> N(l,l) ; N.push_back(2) ; N.push_back(0) ; N.push_back(0) ; // N={1,2,0,0} 

22 i .push_back(SubsystemsInteraction(ptm,N) ) ; 

23 Composite c(i); 

24 

25 // Initial condition 

26 HS_Vector Psi=WavePacket(p,xO,kO,xsig)*Coherent(mPlus, alpha) ^Coherent (mMinus, beta) ; 

27 // Instantiate Trajectory 

28 Trajectory t (Psi , c , seed) ; 

29 // Run Trajectory 

30 RunTrajectory (t ,T) ; 
31 

32 } 

Figure 1: Full driver for one particle in a ring cavity sustaining two travelling- wave modes with opposite 
wave vectors, the ~K mode being pumped. The definition of parameters (omrec, fin, etc.) has been 
omitted for the sake of compactness. 



// Instantiate Frees 

1 LossyMode m(DeltaC, kappa, cutoff) ; // FreeO 

2 PumpedMovingParticle p(omrec , etaef f , f in,K , Sin) ; // Freel, Free2 only one instant! 

3 

4 // Instantiate Interactions 

5 ParticleOrthogonalToCavity pc(m,p,U0); // only one instant! 

6 IdenticalParticles id(p,2,Psileft,Psiright) ; 
7 

8 // Instantiate Composite 

9 std: : vector<SubsystemsInteraction> i (1 , SubsystemsInteraction(pc ,0 , 1) ) ; 

10 i . push_back(SubsystemsInteraction(pc ,0 ,2) ) ; i .push_back(SubsystemsInteraction(id, 1,2)) ; 

11 Composite c(i); 
12 

13 // Initial condition 

14 HS_Vector Psi=Coherent(m,alpha)*TwoParticleState(p,SuperFluid) ; 

Figure 2: The essential part of the driver for two identical pumped particles moving orthogonal to the axis 
of a cavity sustaining one single sinusoidal mode — or otherwise, two identical particles moving in a one 
dimensional optical lattice with the cavity aligned orthogonally to the lattice. 



# MCWFS Driver Parameters: 

# seed=1001 

# eps=le-05 

# dplimit=0.1 

# Displaying in every 10 timestep 

# Composite Dissipative System of Dimension 768 

# Subsystem Nr. 

# Moving Particle: 

# omrec=l 

# Spatial Degree of Freedom finesse=6 



# Subsystem Nr . 1 

# Lossy Mode: 

# Z=(1,0) kappa=l N=3 



# Subsystem Nr . 2 

# Lossy Mode: 

# Z=(1,0) kappa=l N=4 

# eta=(0.3,-0.07) 

# Field from pump: (0.3,-0.07) 

ff 1 <-> Interaction 

ff Particle Moving along Cavity 

# Particle-Cavity Interaction Unot=-l K=l (etaeff =-1) . Mode function type: Plus 

ff 2 <-> Interaction 

# Particle Moving along Cavity 

# Particle-Cavity Interaction Unot=-l K=l (etaeff =-1) . Mode function type: Minus 

# 1 <-> 2 <-> Interaction 

# Particle Two Modes 
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Figure 3: Typical output of the ring-cavity driver of Fig. [H The first two columns are time and time step, 
respectively, then, separated by tab characters, the data stemming from the different subsystems follows: 
columns 3-6 contain the data from subsystem Nr. MovingParticle, columns 7-10 and 11-14 from the two 
cavity modes. The interaction elements make no output in this example. 
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Figure 4: Class inheritance hierarchy starting from the almost purely abstract interface OpenSystem, the 
interface that all simulated systems has to provide for our Trajectory driver. At the bottom of the hierarchy 
we have provided an example set of Elements taken from CQED with moving particles. These may serve as 
building blocks for Composites. The colour code: magenta-framed classes are abstract classes, black-framed 
ones are concrete types; arrows denote class inheritance; in each class the most important functions are 
displayed — purely virtual ones in red, virtual ones in magenta and concrete ones in black; the functions 
displayed in the salmon stripes belong to the private part of the class while the blue and white stripes refer 
to the protected and public part, respectively. The displayed functions are partly documented in the text, 
and partly in the source code. 
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Figure 5: Alternative design featuring a completely separate branch for systems using interaction picture. 
Red-framed classes are virtual bases, and red arrows denote virtual inheritance. 
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Figure 6: The design actually used in the framework. The advantage over the first design is that the 
fundamental functions H, U, J, and Display are declared as pure virtual, and therefore it is very clear 
which class has implemented which function. Still, it does not use virtual bases as the second design. 
The function of Element has ceased to exist so this class is omitted, we have instead a set of classes 
ElementHamiltonian etc. Composite then deals separately with Frees and Interactions. Logically, the 
root of the hierarchy is not called DpenSystem anymore, since the jump function is declared outside this 
class, but merely QuantumSystem. As an example we have plotted LossyMode and PumpedLossyMode to 
show how concrete elements fit into this hierarchy. Note that eg H cannot even be called for LossyMode, 
only for PumpedLossyMode since the first is not derived from the Hamiltonian class. 
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Figure 7: Massive pumped particle moving quickly in a direction orthogonal to the axis of a cavity, (a) 
Expectation value of the atom's position. Each time the atom goes out of the quantisation volume at x = tt, 
it comes back in at x = —ir due to periodic boundary condition, (b) Spread of the atomic wave packet, 
(c) & (d) Real and imaginary part of the scattered field in the cavity, the green lines corresponding to the 
estimation (fT3l). 



State vector viewed via CPA_Slices CPA_View 

1 ... ... 22 23 

FresO 00 00 000 00 000 00 • ■■■■■■ ■ strideS={8) 

Free l 000000000000000000000000 III II strideS={2} 

Free 2 0000 ;0 0000 0000000000 I strideS={l} 

000000000000000000000000 I st rideS= { 8 , 2 } 

0000 000000 00000000 I strideS= { 8, 1 } 

000000000000000000000000 I st rideS= { 2 , 1 } 

000000000000000000000000 I st rideS= { 8 , 2 , 1 } 



Figure 8: The state vector of a system consisting of three subsystems with dimensions 3, 4 and 2, covered 
by different sets of CPA_Slices corresponding to the free subsystems and the interactions between the 
subsystems. One slice is the set of indices displayed in the same colour. A CPA_View is essentially an array 
of slices with the modification that the strides arc stored only once. 



